{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7dfde4cc-a760-4a65-a7ac-50edb273ab74",
   "metadata": {},
   "source": [
    "验证姚端正，数学物理方法，P.282 \n",
    "date: 2024/05/11\n",
    "author: 赵振华"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "id": "b36deca8-56f3-4e1e-92c7-734e1b7f5779",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "import numpy as np\n",
    "from sympy.abc import x, b, l,k\n",
    "c1,c2,c3, c4= symbols('c1 c2 c3 c4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "a3f9064f-dc2e-4e82-81d0-d4f2bc63f9ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = c1*(1-x**2)+c2*(1-x**2)**2+c3*(1-x**2)**3+c4*(1-x**2)**4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "id": "beb8887a-1dc0-4fa4-9371-11995cdb30e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{k \\left(\\frac{c_{1}^{2}}{6} + \\frac{c_{1} c_{2}}{4} + \\frac{c_{1} c_{3}}{5} + \\frac{c_{1} c_{4}}{6} + \\frac{c_{2}^{2}}{10} + \\frac{c_{2} c_{3}}{6} + \\frac{c_{2} c_{4}}{7} + \\frac{c_{3}^{2}}{14} + \\frac{c_{3} c_{4}}{8} + \\frac{c_{4}^{2}}{18}\\right)}{l}$"
      ],
      "text/plain": [
       "k*(c1**2/6 + c1*c2/4 + c1*c3/5 + c1*c4/6 + c2**2/10 + c2*c3/6 + c2*c4/7 + c3**2/14 + c3*c4/8 + c4**2/18)/l"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I = k*integrate (x*y**2,(x,0,1))/l\n",
    "I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "id": "1398bbf4-6014-4fac-9d89-8a4d3b08fa5b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle c_{1}^{2} + \\frac{4 c_{1} c_{2}}{3} + c_{1} c_{3} + \\frac{4 c_{1} c_{4}}{5} + \\frac{2 c_{2}^{2}}{3} + \\frac{6 c_{2} c_{3}}{5} + \\frac{16 c_{2} c_{4}}{15} + \\frac{3 c_{3}^{2}}{5} + \\frac{8 c_{3} c_{4}}{7} + \\frac{4 c_{4}^{2}}{7}$"
      ],
      "text/plain": [
       "c1**2 + 4*c1*c2/3 + c1*c3 + 4*c1*c4/5 + 2*c2**2/3 + 6*c2*c3/5 + 16*c2*c4/15 + 3*c3**2/5 + 8*c3*c4/7 + 4*c4**2/7"
      ]
     },
     "execution_count": 135,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J = integrate (x*diff(y,x)**2,(x,0,1))\n",
    "J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c8b70639-4eab-4255-ad19-5b23677d9079",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 c_{1} + \\frac{4 c_{2}}{3} + c_{3} - k \\left(\\frac{c_{1}}{3} + \\frac{c_{2}}{4} + \\frac{c_{3}}{5}\\right)$"
      ],
      "text/plain": [
       "2*c1 + 4*c2/3 + c3 - k*(c1/3 + c2/4 + c3/5)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Jdc1=diff(J - l*I,c1)\n",
    "Jdc1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "43790510-206b-4228-94b8-071277902fd6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle c_{1} \\cdot \\left(2 - \\frac{k}{3}\\right) + c_{2} \\cdot \\left(\\frac{4}{3} - \\frac{k}{4}\\right) + c_{3} \\cdot \\left(1 - \\frac{k}{5}\\right)$"
      ],
      "text/plain": [
       "c1*(2 - k/3) + c2*(4/3 - k/4) + c3*(1 - k/5)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "collect(expand(Jdc1), (c1,c2,c3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "id": "5391dd2d-08d0-4add-9ac8-ace8d3422123",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0, 0, 0, 0],\n",
       "[0, 0, 0, 0],\n",
       "[0, 0, 0, 0],\n",
       "[0, 0, 0, 0]])"
      ]
     },
     "execution_count": 136,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "matrix_c=zeros(4)\n",
    "matrix_c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "id": "95ca81ba-d663-46b7-b575-5fea769a2eca",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = [c1,c2,c3,c4]\n",
    "\n",
    "for i in [0,1,2,3]:\n",
    "    for j in [0,1,2,3]:\n",
    "        matrix_c[i,j]= expand(diff(J - l*I,c[i])).coeff(c[j])\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "id": "b1692b1e-ff78-4669-a6c1-21b9bd51164c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 - \\frac{k}{3} & \\frac{4}{3} - \\frac{k}{4} & 1 - \\frac{k}{5} & \\frac{4}{5} - \\frac{k}{6}\\\\\\frac{4}{3} - \\frac{k}{4} & \\frac{4}{3} - \\frac{k}{5} & \\frac{6}{5} - \\frac{k}{6} & \\frac{16}{15} - \\frac{k}{7}\\\\1 - \\frac{k}{5} & \\frac{6}{5} - \\frac{k}{6} & \\frac{6}{5} - \\frac{k}{7} & \\frac{8}{7} - \\frac{k}{8}\\\\\\frac{4}{5} - \\frac{k}{6} & \\frac{16}{15} - \\frac{k}{7} & \\frac{8}{7} - \\frac{k}{8} & \\frac{8}{7} - \\frac{k}{9}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[  2 - k/3,   4/3 - k/4,   1 - k/5,   4/5 - k/6],\n",
       "[4/3 - k/4,   4/3 - k/5, 6/5 - k/6, 16/15 - k/7],\n",
       "[  1 - k/5,   6/5 - k/6, 6/5 - k/7,   8/7 - k/8],\n",
       "[4/5 - k/6, 16/15 - k/7, 8/7 - k/8,   8/7 - k/9]])"
      ]
     },
     "execution_count": 138,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "matrix_c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "id": "af96e68f-2a7a-408b-9a64-2fd1573638b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{k^{4}}{10668672000} - \\frac{k^{3}}{27783000} + \\frac{k^{2}}{315000} - \\frac{64 k}{826875} + \\frac{32}{91875}$"
      ],
      "text/plain": [
       "k**4/10668672000 - k**3/27783000 + k**2/315000 - 64*k/826875 + 32/91875"
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "matrix_c.det()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "id": "303fa4c3-fc93-4c05-8d64-a6d95dabb0c9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{5.78319 + 2.3647 \\cdot 10^{-19} i, 30.4747 - 3.25577 \\cdot 10^{-19} i, 78.2312 + 8.08299 \\cdot 10^{-20} i, 269.511 + 8.27757 \\cdot 10^{-21} i\\right\\}$"
      ],
      "text/plain": [
       "{5.78319 + 2.3647e-19*I, 30.4747 - 3.25577e-19*I, 78.2312 + 8.08299e-20*I, 269.511 + 8.27757e-21*I}"
      ]
     },
     "execution_count": 140,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_k = solveset(matrix_c.det(),k)\n",
    "sol_k=sol_k.evalf(6)\n",
    "sol_k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "id": "787b279e-f969-4d9a-9f14-83e6f8e2953a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[  5.78318596],\n",
       "        [ 30.47469711],\n",
       "        [ 78.23121643],\n",
       "        [269.51092529]])"
      ]
     },
     "execution_count": 141,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mat(re(Matrix(list(sol_k))).tolist()).astype(np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baa5998f-3df4-4db3-b3e0-7a7d8b6a5054",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
